function GSS = gs_dataload(filena)
%gs_dataload - GridStrain data loading from an ASCII input file 
%
%       GSS = gs_dataload(filena)
%
% The ASCII file from which the data are read usually contains 10 columns
% (the column number can be 7, see below) in the order
%
%     ns e n ve vn eve evn a b theta
%
% where
% 
%   ns    : station numeric identifier (e.g. 1, 2, ...);
%   e     : east coordinate;
%   n     : north coordinate;
%   ve    : east velocity (or displacement);
%   vn    : north velocity (or displacement);
%   eve   : east velocity (or displacement) error (rms);
%   evn   : north velocity (or displacement) error (rms);
%   a     : length of major half-axis of error ellipse;
%   b     : lenght of minor half-axis of error ellipse;
%   theta : azimuth of major axis of error ellipse,
% 
% or
%
%        statName e n ve vn eve evn a b theta,
%
% where statName(k) is the 4-character string of the k-th station name.
% The function automatically discriminates betweer ns and statName.
%
% The filename of the input ASCII file is filena. If filena is undefined or 
% empty, it can be interactively chosen. The file 'tutorialData2.txt' is
% available for tutorial purposes. 
% The velocity (or displacement) vectors are drawn together with the
% corresponding error ellipses, on the basis of a suitable scaling factor 
% automatically calculated (as usual, 2-sigma error ellipses are shown). 
% 
% If the ASCII file contains 7 columns only, the order is (numeric IDs)
%
%       ns e n ve vn eve evn,
%
% or (IDs represented by strings with station names)
%
%       statName e n ve vn eve evn,
%
% where the variables are as above. In this case, each error ellipse has
% the axes parallel to east and north ones, and the length of half-axes 
% are simply eve and evn. 
%
% Coordinates are in m and velocities, as well as the corresponding errors,
% in mm/y, according to the fact that the velocities are typically obtained
% by GNSS data processing. 
% The function gs_dataload assumes that a mean value of station velocities 
% exceeding 1 surely corresponds to velocities espressed in mm/y. If there
% is the suspect that input velocities are in m/y (e.g., absolute values 
% strongly lower than 1), a dialog box allows the confirmation or exclusion 
% of the fact that the velocities and related data are in m/y instead of 
% mm/y.
%
% After a gn_dataload session, the output struct variable GSS has these 
% fields:
%
%       StationID       (empty if statName is in the input file)
%       StationName     (empty if ns is in the input file)
%       East
%       North
%       EastVelocity
%       NorthVelocity
%       EastVelocityError
%       NorthVelocityError
%       ErrorEllipsea
%       ErrorEllipseb
%       ErrorEllipsetheta
%       RefImageStruct
%
% Note that the output velocities and related uncertainties are always
% expressed in m/y to allow their direct use in strain computation, 
% regardless to the unit used (m/y or mm/y) in the input ASCII file. 
%   
%           GSS = gs_dataload(filena)

% G. Teza, 2005, 2008, 2022.

Options = GridStrainOptions;

fsz = 14;       % font size 
fszr = 10;      % reduced font size
LW = 1.5;

if (nargin < 1) || ~ischar(filena)
    filena = [];              
end

format long e

nofi = 0;
while nofi == 0
    if isempty(filena)
        [filen, pathn] = uigetfile(...
            {'*.txt','ASCII (*.txt) files';...
            '*.dat','ASCII (*.dat) files';...
            '*.*','ASCII (*.*) files'},...
            'INPUT ASCII FILENAME');
        if isequal(filen,0) || isequal(pathn,0)
            filena = [];    % cancel was pressed by user
        else
            filena = fullfile(pathn,filen);
        end
    end
    if ~isempty(filena)
        try
            % case numerical station ID
            A = load(filena);
            statName = [];
            nvar = size(A,2);
        catch
            % case station name
            fid = fopen(filena);
            CA = textscan(fid,'%s%f%f%f%f%f%f%f%f%f');
            fclose(fid);
            CAS = CA(1); 
            statName = CAS{:};
            ca8 = CA{8};    % test on 8-th column
            nrs = numel(ca8);
            if sum(isnan(ca8)) == nrs
                CA(8:end) = []; % in this case, no error ellipses are provided
            end
            ANC = CA(2:end);
            ncn = size(ANC,2);
            AN = nan(nrs,ncn);
            for k = 1:ncn 
                AN(:,k) = ANC{k};
            end
            A = [(1:nrs)' AN];
            nvar = ncn+1;
        end
        Ina = sum(isnan(A),2) > 0;  % Suppression of rows having at least a NaN element 
        A(Ina,:) = [];
        if ~isempty(statName)
            statName(Ina) = []; %#ok<AGROW>
        end
    end
    if isempty(A)
        strFileno = 'INVALID OR EMPTY FILE'; 
    elseif size(A,1) < 2
        strFileno = 'NOT ENOUGH EXPERIMENTAL POINTS';
    elseif nvar < 7
        strFileno = 'THE COLUMNS MUST BE AT LEAST 7';
    else
        strFileno = '';
        nofi = 1;
    end
    if ~isempty(strFileno)
        [optExit,v] = listdlg(...
            'PromptString',sprintf('INVALID CHOICE: %s',strFileno),...
            'SelectionMode','single',...
            'ListString',{'ANOTHER CHOICE','EXIT'},...
            'ListSize',[300 200]);
        if ~v || (optExit == 2)
            nofi = 1;
            A = [];
            GSS = emptyout;     % ancillary function, see below
        end
    end
end

if isempty(A)
    return
end

if isempty(statName)
    ns = A(:,1); 
else
    ns = [];
end
e  = A(:,2); 
n  = A(:,3);
lv = length(e);

%-------------------------------------------------------------------------
[ve,nido] = mm2m(A(:,4),2,5,'East');  % Ancillary function, see below
vn  = mm2m(A(:,5),nido);
eve = mm2m(A(:,6),nido);
evn = mm2m(A(:,7),nido);
scalef = 0.2*max(max(e)-min(e),max(n)-min(n))/max(max(abs(ve),abs(vn)));
%-------------------------------------------------------------------------

if nvar < 10
    a = 2*eve*scalef;
    b = 2*evn*scalef;
    theta = zeros(lv,1);
else
    a = 2*mm2m(A(:,8),nido)*scalef;
    b = 2*mm2m(A(:,9),nido)*scalef;
    theta = pi/180*A(:,10);
end

%% Possible choice of an error threshold 

sqab = sqrt(eve.^2+evn.^2);   % Square sum of a and b
hki = figure; 
set(hki,'Units','normalized','OuterPosition',[0.1 0.1 0.8 0.8]);
mab = histnorm(sqab);
xlabel('Uncertainty distribution','FontSize',fsz); 

lbi1 = uicontrol(hki,'Style','listbox',...
    'string', {'DEFINE AN ERROR THRESHOLD','SKIP'},...
    'FontSize',fsz,...
    'Units','normalized','Position', [0.05 0.05 0.30 0.15],...
    'Callback', 'uiresume(gcbf)');
    uiwait(hki);
    vlbi1 = lbi1.Value;
    delete(lbi1);

if vlbi1 == 1
    prompt = {'Threshold (mm): '}; 
    defans = {num2str(mab)}; 
    fields = {'thre'};
    nocont = 0;
    while nocont == 0
        acquis = inputdlg(prompt,'THRESHOLD',1,defans);
        if ~isempty(acquis)
            acq  = cell2struct(acquis,fields);
            thre = str2double(acq.thre);
            if ~isempty(thre) && (thre > 0) 
                lbi2 = uicontrol(hki,...
                    'Style', 'listbox',...
                    'string',...
                        {sprintf('ACCEPT THRESHOLD %f mm',thre),...
                        'SELECT A DIFFERENT VALUE',...
                        'SKIP THRESHOLD CHOICE'},...
                    'FontSize',fsz,...
                    'Units','normalized',...
                    'Position', [0.05 0.05 0.30 0.15],...
                    'Callback', 'uiresume(gcbf)');
                    uiwait(hki);
                    vlbi2 = lbi2.Value;
                    delete(lbi2);
                    menacc = vlbi2;
                    iok = 1;
            else
                iok = 0;
            end
        else
            iok = 0;
        end
        if ~iok
            lbi3 = uicontrol(hki,...
                'Style', 'listbox',...
                'string',...
                    {'INVALID INPUT - NEW SELECTION','SKIP THRESHOLD CHOICE'},...
                'FontSize',fsz,...
                'Units','normalized',...
                'Position', [0.05 0.05 0.30 0.15],...
                'Callback', 'uiresume(gcbf)');
                uiwait(hki);
                vlbi3 = lbi3.Value;
                delete(lbi3);
                menacc = vlbi3+1;
        end
        if menacc == 1
            nocont = 1;
        elseif menacc == 3
            nocont = 1; 
            thre = [];
        end
    end
    if ~isempty(thre)
        Ift = sqab < thre;
        eve(Ift) = thre;
        evn(Ift) = thre; 
        a(Ift)   = thre*scalef; 
        b(Ift)   = thre*scalef;
        theta(Ift) = 0;
    end
end
close(hki);

GSS.StationID = ns;
GSS.StationName = statName;
GSS.East = e;
GSS.North = n;
GSS.EastVelocity = ve;
GSS.NorthVelocity = vn;
GSS.EastVelocityError = eve;
GSS.NorthVelocityError = evn;
GSS.ErrorEllipsea = a;
GSS.ErrorEllipseb = b;
GSS.ErrorEllipsetheta = theta;

%% DATA VISUALIZATION

%---- for the "arrow scale" visualization 
esca = e;  
nsca = n;
vesca = ve; 
vnsca = vn;
% A new point is generated in order to represent a velocity vector along E axis
std_esca = std(esca-mean(esca));
std_nsca = std(nsca-mean(nsca));
esca(lv+1)  = min(esca)-std_esca/2;
nsca(lv+1)  = min(nsca)-std_nsca/2;
vesca(lv+1) = (max(abs(ve))+max(abs(vn)))*0.75;
valnum = vesca(end);
if nido == 1
    valnum = 1000*valnum; 
end
valnapp = ApprInt(valnum);
val = num2str(valnapp);
if nido == 1
    valnapp = valnapp/1000; 
end
vesca(end) = valnapp;  % in this way, an optimal vector is drawn 
vnsca(lv+1) = 0.00;

% visualization:
figki = figure;
set(figki,'Units','normalized','OuterPosition',[0 0 1 1]); 
hold on; 
plot(e,n,'xr','LineWidth',LW);
quiver(esca,nsca,vesca*scalef,vnsca*scalef,0,...
    'b','LineWidth',LW); % vector visualization 

for k = 1:lv
    [xell,yell] = ellipse(2*a(k),2*b(k),...
        e(k)+vesca(k)*scalef,n(k)+vnsca(k)*scalef,theta(k),100);
    plot(xell,yell,'r');
    if Options.VisualizePointIDs
        se = sign(vesca(k)); 
        sn = sign(vnsca(k));
        if isempty(statName)
            strStatk = num2str(ns(k));
        else
            strStatk = statName{k};
        end
        ht = text(e(k)-se*std_nsca/15,n(k)-sn*std_nsca/15,strStatk);
        set(ht,'Color','k','FontSize',fszr);
    end
end
htg = text(esca(end),nsca(end)+std_nsca/10,[val ' mm/y']);
set(htg,'FontSize',fszr);
title('VELOCITIES AND ERROR ELLIPSES','FontSize',fsz);
set(gca,'FontSize',fszr);
xlabel('EAST (m)','FontSize',fsz); 
ylabel('NORTH (m)','FontSize',fsz);
axis equal;
ALA = axis; 

%% POSSIBLE BACKGROUND IMAGE
lbgim1 = uicontrol(figki,'Style', 'listbox',...
    'string', {'USE A BACKGROUND IMAGE','SKIP (NO BACKGROUND IMAGE)'},...
    'FontSize',fsz,...
    'Units','normalized','Position', [0.05 0.05 0.30 0.15],...
    'Callback', 'uiresume(gcbf)');
    uiwait(figki);
    vlbgim1 = lbgim1.Value;
    delete(lbgim1);
if vlbgim1 == 1
    niter = 0;
    while niter == 0
        [fileim,pathim] = uigetfile(...
            {'*.tif;*.tiff;*.jpg;*.jpeg',...
            'Image files (*.tif,*.tiff,*.jpg,*.jpeg)'},...
            'INPUT IMAGE FILE OR MATLAB MANAGING FILE');
        if ~isequal(fileim,0) && ~isequal(pathim,0)
            fileIm = fullfile(pathim,fileim);
            [~,fim,~] = fileparts(fileIm);
            fileVert = fullfile(pathim,[fim '.mat']);
            if exist(fileVert,'file') 
                niter = 1;
                IOk = 1;
            else
                IOk = 2;
                strInv = 'INVALID MATLAB VERTICES FILE';
            end
        else
            IOk = 3;
            strInv = 'INVALID IMAGE';
        end
        if IOk > 1 
            lbgim2 = uicontrol(figki,'Style','listbox',...
                'string',{sprintf('%s - NEW CHOICE',strInv),...
                    'SKIP (NO BACKGROUND IMAGE)'},...
                'FontSize',fsz,...
                'Units','normalized','Position', [0.05 0.05 0.30 0.15],...
                'Callback', 'uiresume(gcbf)');
            uiwait(figki);
            vlbgim2 = lbgim2.Value;
            delete(lbgim2);
            if vlbgim2 == 2
                fileIm = [];
                niter = 1; 
            end
        end
    end
    if ~isempty(fileIm)
        IMS = struct('Image',fileIm,'Vertices',fileVert);
    else
        IMS = [];
    end
else
    IMS = [];
end
gs_ReadBGImage(IMS,ALA);
hold off;
axis equal; 

GSS.RefImageStruct = IMS;

%% ANCILLARY FUNCTIONS

function GSS = emptyout
GSS = struct('StationID',[],'StationName',[],'East',[],'North',[],...
    'EastVelocity',[],'NorthVelocity',[],...
    'EastVelocityError',[],'NorthVelocityError',[],...
    'ErrorEllipsea',[],'ErrorEllipseb',[],'ErrorEllipsetheta',[],'IMS',[]);

function [vdd,nido] = mm2m(vd,nidi,tresh,strX)
% If nidi is 0, vdd = vd and nido = 0.
% If nidi is 1, vdd = vd/1000 and nido = 1.
% If nidi is 2 and max(|vd|) > tresh, vdd = vd and nido = 1 (if thresh is
% undefined or empty, it is thresh = Inf).
% If nidi is 2 and max(|vd|) <= tresh, a warning message is emitted because 
% it seems that velocities are expressed in m/y (not in mm/y, as expected). 
% In this case, the user can iteractively confirm on a velocity histogram
% that the velocities are expressed in mm/y or say that they are in m/y. 
% If they are in mm/y it is vdd = vd and nido = 0. If they are in m/y it is
% vdd = vdd/1000 and nido = 1, to have % always velocities in mm/y.
% If nidi is 2 and mean(x) > tresh, vdd = vd/1000 and nido = 1.
% The optional input arguments strX applies only if the histogram is shown.
% This string strX is placed before 'velocity distribution' in the xlabel.
% If undefined or empty, it is strX = ''. 
if (nargin < 4) || isempty(strX)
    strX = '';
end
if (nargin < 3) || isempty(tresh)
    tresh = Inf; 
end
if nidi == 0
    vdd = vd; 
    nido = 0; 
end
if nidi == 1
    vdd = vd/1000; 
    nido = 1; 
end    
if nidi == 2
    vema = max(abs(vd));
    if vema >= tresh
        vdd = vd/1000; 
        nido = 1; 
    else
        hmm = figure;
        set(hmm,'Units','normalized','OuterPosition',[0.1 0.1 0.8 0.8]);
        histnorm(vd);
        xlabel(sprintf('%s velocity distribution',strX),'FontSize',14); 
        lbmm = uicontrol(hmm,'Style','listbox',...
        'string',...
            {sprintf('WARNING: MAX VELOCITY IS %f mm/y; IS IT REALLY IN mm/y?',vema),...
            'NO, CONVERSION m -> mm IS NECESSARY'},...
        'FontSize',14,...
        'Units','normalized','Position', [0.05 0.05 0.30 0.15],...
        'Callback', 'uiresume(gcbf)');
        uiwait(hmm);
        vlbmm = lbmm.Value;
        delete(lbmm);
        if vlbmm == 1
            vdd = vd/1000; 
            nido = 1;
        else
            vdd = vd; 
            nido = 0;
        end
        close(hmm);
    end
end